This code runs the analyses in the paper. Below we lay out the process in detail for one dataset, and then provide the code for the other datasets below this. Note that much of this uses functions in the simulations.R file.

set.seed(123)

Load packages and functions

Make sure you work with dispRity 0.4 and TG’s version of Claddis (much faster than GL’s).

# The latest version of devtools and ape
library(devtools)
library(ape)

## The latest version of Claddis (on development on TG branch)
## install_github("TGuillerme/Claddis")
library(Claddis)

## The latest version of dispRity on master branch
## install_github("TGuillerme/dispRity", ref = "master")
library(dispRity) 

## Load functions using source
# Doesn't always work... Rmd problem?
source("../functions/simulations_fun.R")

Loading the data

To run dispRity analyses we need the morphospace, the phylogeny (with a root age and node labels) and the first and last occurrence dates (FADLAD). The tree and FADLAD are in the data/ folder. To construct the morphospace we need to use the distance_matrix files we created in 01-extract-data-for analyses.Rmd and stored in data/processed.

## Choose a slug
slug <- "Beck2014"

## Loading the discrete morphological matrix
## This is just to match the names with the tree
matrix <- ReadMorphNexus(paste0("../data/matrices/", slug, ".nex"))

## Loading the FADLAD
FADLAD <- read.csv(paste0("../data/FADLAD/", slug, ".csv"), row.names = 1, header = TRUE)

## Loading the tree
tree <- read.nexus(paste0("../data/trees/", slug, ".tre"))
## Adding node labels and a root age (max tree age date)
tree <- makeNodeLabel(tree, method = "number", prefix = "n")
tree$root.time <- max(tree.age(tree)[, 1])

## Cleaning the tree to remove things that aren't in the matrix 
cleaned_data <- clean.data(matrix$matrix, tree)
matrix$matrix <- cleaned_data$matrix
tree <- cleaned_data$tree

## Loading the distance matrices
## See "01-extract-data-for analyses.Rmd" for how we made these
load(paste0("../data/processed/distance_matrix.", slug ,".Rda"))

## Ordinating the matrix
morphospace <- cmdscale(matrix_dist, k = nrow(matrix_dist) - 2, add = TRUE)$points

Scenarios

We used the data to explore these different scenarios:

  1. Stratigraphy. This is the traditional method, where all the species within each stratigraphic period are included in the disparity calculation. This often leads to bins of unequal duration.

  2. Equally sized time bins.
    1. Where the duration of the bin is equal to the median duration of the stratigraphic period.
    2. Where the number of bins is equal to the number of stratigraphic periods.
    3. Where the number of bins is either 2, 5, 10, 15, or 20.
  3. Time slicing.
    1. Where the interval between slices is equal to the median duration of the stratigraphic period.
    2. Where the number of slices is equal to the number of stratigraphic periods.
    3. Where the number of slices is either 2, 5, 10, 15, or 20.

Running the analyses and saving the outputs

These analyses use a series of functions in functions\. These essentially wrap the existing dispRity functions so that running lots of models and extracting lots of data is not such a big issue.

TG: for the following pipeline, I’ve reduced the number of rarefaction that highly increased the calculation time. We now only calculate the disparity for the non-rarefied subsamples and the subsamples rarefied to contain only 3 taxa (the lowest level of rarefaction). We can use the full rarefaction if we want to plot rarefaction curves later on.

## Run the function for Beck data
slug <- "BeckLee"

## Ages (rarefied to 3 taxa only)
out1 <- run.all.disparity.wrapper(morphospace, tree, type = "Age", FADLAD, inc.nodes = TRUE,
                                bootstraps = 100, metric = c(sum, variances),
                                metric.name = "sum_var", rarefaction = FALSE)
#readr::write_csv(out1, path = paste0("../outputs/age", "_", slug, ".csv"))

## Epochs (rarefied to 3 taxa only)
out2 <- run.all.disparity.wrapper(morphospace, tree, type = "Epoch", FADLAD, inc.nodes = TRUE,
                        bootstraps = 100, metric = c(sum, variances), metric.name = "sum_var",
                        rarefaction = FALSE)
#readr::write_csv(out2, path = paste0("../outputs/epoch", "_", slug, ".csv"))

Note that you can use the previous output format (the data tables) as before by accessing them through out$results$stratigraphy/duration/number. Writing them as a csv will require some thinking since they don’t have the same number of columns (though it’s still doable I guess). Also, if you want to plot things through ggplot, you can use the values in these data tables!

Examining the output

The results as a table

Disparity per stratigraphic time: unequal time bins and non-equidistant time slices
subsamples n obs bs.median 2.5% 25% 75% 97.5% model stratigraphy bin_type
163.5 - 145 10 33.32281 30.22177 26.63719 29.25598 31.10961 31.97925 equalbins Epoch stratigraphy
145 - 100.5 32 34.05927 33.07039 32.21348 32.81688 33.31503 33.65910 equalbins Epoch stratigraphy
100.5 - 66 82 35.03925 34.61814 34.29401 34.49754 34.72298 34.91502 equalbins Epoch stratigraphy
66 - 56 35 35.75072 34.77715 34.01151 34.61513 34.99665 35.32506 equalbins Epoch stratigraphy
56 - 33.9 35 35.98760 35.02616 34.24015 34.81680 35.21084 35.51950 equalbins Epoch stratigraphy
33.9 - 23.03 3 37.40622 24.97051 0.00000 24.48371 25.35822 37.40622 equalbins Epoch stratigraphy
23.03 - 5.332 1 NA NA NA NA NA NA equalbins Epoch stratigraphy
5.332 - 2.58 0 NA NA NA NA NA NA equalbins Epoch stratigraphy
2.58 - 0 14 36.58173 33.99437 31.26595 33.19363 34.56647 35.49716 equalbins Epoch stratigraphy
172.75 0 NA NA NA NA NA NA acctran Epoch stratigraphy
167.25 3 33.66778 23.38683 0.00000 21.46503 23.38683 33.66778 acctran Epoch stratigraphy
117.75 15 34.93270 32.72896 30.21917 32.08867 33.35239 33.96271 acctran Epoch stratigraphy
71 43 35.67388 34.88899 34.12592 34.69091 35.04140 35.28352 acctran Epoch stratigraphy
67.05 45 35.76368 34.98565 34.34750 34.89697 35.13146 35.44430 acctran Epoch stratigraphy
39.335 23 36.71698 35.09449 33.77583 34.81015 35.44410 35.93741 acctran Epoch stratigraphy
31.879 17 36.73472 34.67493 32.91496 34.07698 35.12814 35.69984 acctran Epoch stratigraphy
6.708 14 36.58173 34.12833 31.31029 33.51234 34.67076 35.74088 acctran Epoch stratigraphy
3.87 14 36.58173 34.07553 31.72960 33.37189 34.58847 35.41521 acctran Epoch stratigraphy
172.75 0 NA NA NA NA NA NA deltran Epoch stratigraphy
167.25 2 31.10128 15.55064 0.00000 0.00000 31.10128 31.10128 deltran Epoch stratigraphy
117.75 12 33.05388 30.63674 28.21955 29.77720 31.21053 31.87566 deltran Epoch stratigraphy
71 29 34.14936 33.05078 32.22101 32.65106 33.26051 33.54114 deltran Epoch stratigraphy
67.05 29 34.16739 32.95663 31.95007 32.64143 33.26758 33.63848 deltran Epoch stratigraphy
39.335 15 33.29550 31.19177 29.28215 30.51424 31.62102 32.41087 deltran Epoch stratigraphy
31.879 11 32.85585 29.96998 27.59435 29.35261 30.59878 31.58435 deltran Epoch stratigraphy
6.708 9 32.49293 29.64869 26.02903 28.29220 30.61162 31.58760 deltran Epoch stratigraphy
3.87 9 32.49293 28.96746 24.67337 27.98939 30.18050 31.57041 deltran Epoch stratigraphy
172.75 0 NA NA NA NA NA NA punctuated Epoch stratigraphy
167.25 2 33.74351 0.00000 0.00000 0.00000 33.74351 33.74351 punctuated Epoch stratigraphy
117.75 14 33.68405 31.39937 28.91869 30.72099 31.98270 32.58289 punctuated Epoch stratigraphy
71 39 34.86233 34.01326 33.31116 33.78864 34.25030 34.56706 punctuated Epoch stratigraphy
67.05 42 34.88953 34.08132 33.41392 33.87639 34.26933 34.61587 punctuated Epoch stratigraphy
39.335 22 34.72135 33.20598 32.18676 32.65332 33.60269 34.25076 punctuated Epoch stratigraphy
31.879 15 34.99709 32.68275 30.40956 32.11415 33.42262 34.12985 punctuated Epoch stratigraphy
6.708 13 34.60840 32.27409 29.12841 31.39463 32.83334 33.75597 punctuated Epoch stratigraphy
3.87 12 34.64232 32.24175 29.23584 31.27131 32.76387 33.88411 punctuated Epoch stratigraphy
172.75 0 NA NA NA NA NA NA gradual Epoch stratigraphy
167.25 3 31.46531 21.09128 0.00000 20.73418 21.10516 31.46531 gradual Epoch stratigraphy
117.75 14 33.88029 31.69960 29.99147 30.95951 32.19942 33.05888 gradual Epoch stratigraphy
71 32 34.64399 33.54515 32.39444 33.28298 33.83946 34.24164 gradual Epoch stratigraphy
67.05 38 34.99828 34.12221 33.49100 33.94587 34.27612 34.69803 gradual Epoch stratigraphy
39.335 19 34.83142 33.14431 31.33004 32.71437 33.48554 34.10599 gradual Epoch stratigraphy
31.879 14 35.57663 33.05639 30.48743 32.37587 33.79526 34.72970 gradual Epoch stratigraphy
6.708 14 36.58173 34.16908 31.41381 33.45433 34.75614 35.54902 gradual Epoch stratigraphy
3.87 14 36.58173 34.15189 32.15575 33.61581 34.65749 35.41226 gradual Epoch stratigraphy
Disparity per average stratigraphic time: equal time bins and equidistant time slices based on the median stratigraphic epoch duration
subsamples n obs bs.median 2.5% 25% 75% 97.5% model stratigraphy bin_type
159.282 - 141.584 7 33.56010 28.82605 22.82095 27.08952 30.63607 32.19139 equalbins Epoch duration
141.584 - 123.886 15 34.26321 31.92322 29.95243 31.25537 32.51200 33.47121 equalbins Epoch duration
123.886 - 106.188 14 34.22967 31.92039 29.67577 31.04999 32.64491 33.22078 equalbins Epoch duration
106.188 - 88.49 33 34.65900 33.75171 32.99633 33.50910 33.89704 34.17113 equalbins Epoch duration
88.49 - 70.792 45 34.89283 34.14289 33.48815 33.93394 34.33884 34.58966 equalbins Epoch duration
70.792 - 53.094 52 35.61019 34.95198 34.31254 34.74795 35.09250 35.34162 equalbins Epoch duration
53.094 - 35.396 23 35.98965 34.52512 33.51641 34.24298 34.88787 35.29425 equalbins Epoch duration
35.396 - 17.698 4 37.73019 31.45205 18.36278 31.06666 31.77800 37.73019 equalbins Epoch duration
17.698 - 0 15 36.64181 34.18229 31.50106 33.61830 34.78979 35.28545 equalbins Epoch duration
168.131 3 33.66778 22.48369 0.00000 21.46503 23.38683 33.66778 acctran Epoch duration
150.433 6 34.30504 28.57096 19.43692 26.51226 30.59620 32.39496 acctran Epoch duration
132.735 14 35.00759 32.63174 30.29309 31.90464 33.29958 34.05483 acctran Epoch duration
115.037 15 34.93270 32.78958 30.95106 32.11680 33.29410 33.87637 acctran Epoch duration
97.339 28 35.34094 34.15773 33.01856 33.87996 34.32740 34.59344 acctran Epoch duration
79.641 36 35.27742 34.29577 33.32058 34.04266 34.51080 34.95222 acctran Epoch duration
61.943 37 35.82991 34.92185 33.94244 34.69257 35.07417 35.43030 acctran Epoch duration
44.245 23 36.53531 35.02771 33.96363 34.76890 35.28234 35.80668 acctran Epoch duration
26.547 15 36.64181 34.42214 32.18929 33.74064 34.84109 35.45842 acctran Epoch duration
168.131 2 31.10128 0.00000 0.00000 0.00000 31.10128 31.10128 deltran Epoch duration
150.433 3 32.06333 21.64125 0.00000 20.73418 21.75122 32.06333 deltran Epoch duration
132.735 10 32.49117 29.37183 26.01487 28.61740 30.21673 31.13435 deltran Epoch duration
115.037 12 33.05388 30.32200 28.34990 29.66091 31.06855 31.89017 deltran Epoch duration
97.339 21 33.54587 31.95197 31.06267 31.68536 32.30157 32.69440 deltran Epoch duration
79.641 26 34.24332 32.96612 31.66601 32.72204 33.25548 33.70782 deltran Epoch duration
61.943 23 33.79279 32.37213 31.13411 31.97450 32.70949 33.23975 deltran Epoch duration
44.245 16 33.21886 31.22745 29.57337 30.77644 31.89136 32.33618 deltran Epoch duration
26.547 10 32.74774 29.83877 26.42110 28.81568 30.49635 31.31324 deltran Epoch duration
168.131 2 31.10128 0.00000 0.00000 0.00000 31.10128 31.10128 punctuated Epoch duration
150.433 6 34.11709 29.41505 20.82254 26.90145 30.59602 32.53259 punctuated Epoch duration
132.735 14 34.31296 31.86428 29.46666 31.20959 32.52878 33.47351 punctuated Epoch duration
115.037 15 33.89808 31.76499 30.21881 31.25125 32.29596 33.02712 punctuated Epoch duration
97.339 24 34.15318 32.79371 31.92913 32.45932 33.06354 33.54609 punctuated Epoch duration
79.641 34 34.63371 33.66572 32.92233 33.38249 33.90457 34.15045 punctuated Epoch duration
61.943 36 35.01892 34.13122 33.35926 33.89556 34.42769 34.77053 punctuated Epoch duration
44.245 20 34.75993 33.08432 31.50550 32.66992 33.44955 34.00155 punctuated Epoch duration
26.547 15 36.04577 33.74115 31.72936 33.16945 34.18529 35.14037 punctuated Epoch duration
168.131 3 31.46531 21.09128 0.00000 20.73418 31.46531 31.46531 gradual Epoch duration
150.433 5 32.80692 26.10963 16.32335 23.38540 29.42816 29.80455 gradual Epoch duration
132.735 10 33.50439 30.17059 27.73548 29.32479 31.34921 32.50129 gradual Epoch duration
115.037 15 34.03641 31.78850 30.00475 31.24314 32.46321 33.15886 gradual Epoch duration
97.339 24 34.73893 33.31136 32.20370 32.99455 33.58992 34.07044 gradual Epoch duration
79.641 31 34.79364 33.57059 32.74380 33.24192 33.83250 34.26474 gradual Epoch duration
61.943 31 34.86539 33.78412 33.15643 33.54985 34.07509 34.48359 gradual Epoch duration
44.245 18 34.66605 32.86028 31.07256 32.18227 33.36466 34.20057 gradual Epoch duration
26.547 13 36.24387 33.55574 30.78633 32.67200 34.42692 35.16849 gradual Epoch duration
Disparity per average stratigraphic time: equal time bins and equidistant time slices based on the number of stratigraphic epochs
subsamples n obs bs.median 2.5% 25% 75% 97.5% model stratigraphy bin_type
168.35028 - 149.644693333333 7 33.09887 28.59308 22.17681 27.35666 30.31330 32.16348 equalbins Epoch number
149.644693333333 - 130.939106666667 10 32.82549 29.76916 27.28616 28.69485 30.45846 31.89242 equalbins Epoch number
130.939106666667 - 112.23352 16 34.64608 32.67191 30.80198 32.03644 33.15730 33.75325 equalbins Epoch number
112.23352 - 93.5279333333333 16 34.19051 32.03013 30.34755 31.59049 32.61587 33.41532 equalbins Epoch number
93.5279333333333 - 74.8223466666667 53 34.85964 34.23060 33.70008 34.08296 34.37971 34.61608 equalbins Epoch number
74.8223466666667 - 56.11676 54 35.43996 34.78783 34.35541 34.60311 34.91857 35.17115 equalbins Epoch number
56.11676 - 37.4111733333333 28 35.80694 34.59685 33.70423 34.26161 34.85875 35.24247 equalbins Epoch number
37.4111733333333 - 18.7055866666667 9 37.03233 33.12955 28.23091 31.68900 34.75597 35.96668 equalbins Epoch number
18.7055866666667 - 0 15 36.64181 34.49749 32.57691 33.97456 34.87068 35.63350 equalbins Epoch number
168.35028 3 33.66778 22.48369 0.00000 21.46503 23.38683 33.66778 acctran Epoch number
147.306495 8 34.58133 30.72673 25.95108 29.31449 31.81122 33.35937 acctran Epoch number
126.26271 15 34.91262 32.52210 30.97702 31.84374 33.22136 34.20434 acctran Epoch number
105.218925 22 35.04096 33.50811 32.09212 33.26416 33.84137 34.20487 acctran Epoch number
84.17514 32 35.25581 34.20205 33.39683 33.95360 34.42747 34.79412 acctran Epoch number
63.131355 40 35.86652 35.00329 34.43719 34.77424 35.16076 35.43566 acctran Epoch number
42.08757 23 36.71698 35.14141 34.18355 34.78785 35.57226 36.00312 acctran Epoch number
21.043785 15 36.64181 34.28187 32.08969 33.59899 34.74871 35.46550 acctran Epoch number
0 14 36.58173 34.17409 31.74489 33.49142 34.75215 35.41788 acctran Epoch number
168.35028 2 31.10128 0.00000 0.00000 0.00000 31.10128 31.10128 deltran Epoch number
147.306495 4 32.16469 26.87273 15.81309 21.09545 26.90892 32.16469 deltran Epoch number
126.26271 11 32.46694 29.70603 26.23159 28.58966 30.53560 31.22416 deltran Epoch number
105.218925 17 33.18290 31.14234 29.65922 30.48352 31.62709 32.27118 deltran Epoch number
84.17514 20 33.86498 32.32130 30.98335 31.84598 32.55487 33.07225 deltran Epoch number
63.131355 24 33.70318 32.29995 31.40135 32.04073 32.64849 33.22530 deltran Epoch number
42.08757 15 33.29550 31.21691 28.56696 30.62232 31.65407 32.26817 deltran Epoch number
21.043785 10 32.74774 29.73532 25.24299 28.66892 30.67990 31.48521 deltran Epoch number
0 14 36.58173 33.95306 31.82532 33.37233 34.46715 35.36812 deltran Epoch number
168.35028 3 33.48689 22.49567 0.00000 21.09128 23.38683 33.48689 punctuated Epoch number
147.306495 6 32.76166 27.98110 19.07316 25.39128 28.94764 30.79350 punctuated Epoch number
126.26271 12 33.87153 31.35454 28.99568 30.53696 31.90175 32.84359 punctuated Epoch number
105.218925 20 33.89993 32.37231 31.09740 31.91930 32.77386 33.26793 punctuated Epoch number
84.17514 32 34.74205 33.78045 33.00744 33.45730 33.97628 34.36286 punctuated Epoch number
63.131355 39 35.09481 34.25245 33.61823 34.01326 34.48994 34.75350 punctuated Epoch number
42.08757 22 35.59579 33.90372 32.82590 33.62971 34.43416 34.86670 punctuated Epoch number
21.043785 14 34.98620 32.66405 30.37380 31.91377 33.15661 33.89102 punctuated Epoch number
0 14 36.58173 33.99161 31.77111 33.34211 34.57252 35.36994 punctuated Epoch number
168.35028 2 31.10128 31.10128 0.00000 0.00000 31.10128 31.10128 gradual Epoch number
147.306495 5 32.80692 25.93453 19.26351 22.79376 28.89722 30.07253 gradual Epoch number
126.26271 12 33.24868 30.67097 27.71580 29.83864 31.23781 31.98489 gradual Epoch number
105.218925 19 34.07954 32.45076 30.75566 31.89571 32.74788 33.22262 gradual Epoch number
84.17514 25 34.81277 33.55705 32.13463 33.19140 33.85160 34.30737 gradual Epoch number
63.131355 33 35.00839 33.97028 32.88982 33.71818 34.21615 34.57940 gradual Epoch number
42.08757 17 34.57695 32.60078 30.73517 31.87518 33.11267 33.85894 gradual Epoch number
21.043785 14 36.49673 33.95655 32.19033 33.42689 34.45953 35.13698 gradual Epoch number
0 14 36.58173 34.25095 32.40351 33.66960 34.79317 35.35109 gradual Epoch number

This is not really digest and we should find a better way to present that. Anyway, these boring tables are always good for the supplementaries.

Plotting the results dispRity style

Here is an example of the results using plot.dispRity. The upper panel is with the CIs and the lower one is the same without the CIs.

The results for the Age type

op <- par(bty = "n", mfrow = c(3,2))
## Principal colors
colors <- palette()[1:5][-1] #TG: that can be made fancier
colors[5] <- "black"
## Color combination (for the CI)
colors_CI <- c("pink", "lightgreen", "lightblue", "lightcyan", "lightgray")
## y limits
ylim <- range(unlist(lapply(out1$results, function(X) return(c(min(X[,5], na.rm = TRUE), max(X[,5],
  na.rm = TRUE))))))
## Zooming
ylim_zoom <- c(20, 35)

## Plotting the different methods (with/without CIs)
plot.results.dispRity(out1$object, method = 1, colors = colors, colors_CI = colors_CI,
  main = "Stratigraphy", ylim = ylim, CI = TRUE, legend = TRUE)
plot.results.dispRity(out1$object, method = 1, colors = colors, colors_CI = colors_CI,
  main = "Stratigraphy (no CI)", ylim = ylim_zoom, CI = FALSE, time.bins.line = TRUE)

plot.results.dispRity(out1$object, method = 2, colors = colors, colors_CI = colors_CI,
  main = "Duration", ylim = ylim, CI = TRUE)
plot.results.dispRity(out1$object, method = 2, colors = colors, colors_CI = colors_CI,
  main = "Duration (no CI)", ylim = ylim_zoom, CI = FALSE, time.bins.line = TRUE)

plot.results.dispRity(out1$object, method = 3, colors = colors, colors_CI = colors_CI,
  main = "Number", ylim = ylim, CI = TRUE)
plot.results.dispRity(out1$object, method = 3, colors = colors, colors_CI = colors_CI,
  main = "Number (no CI)", ylim = ylim_zoom, CI = FALSE, time.bins.line = TRUE)

The results for the Epoch type

op <- par(bty = "n", mfrow = c(3,2))
## Principal colors
colors <- palette()[1:5][-1] #TG: that can be made fancier
colors[5] <- "black"
## Color combination (for the CI)
colors_CI <- c("pink", "lightgreen", "lightblue", "lightcyan", "lightgray")
## y limits
ylim <- range(unlist(lapply(out2$results, function(X) return(c(min(X[,5], na.rm = TRUE),
  max(X[,5], na.rm = TRUE))))))

## Zooming
ylim_zoom <- c(20, 35)

## Plotting the different methods (with/without CIS)
plot.results.dispRity(out2$object, method = 1, colors = colors, colors_CI = colors_CI,
  main = "Stratigraphy", ylim = ylim, CI = TRUE, legend = TRUE)
plot.results.dispRity(out2$object, method = 1, colors = colors, colors_CI = colors_CI,
  main = "Stratigraphy (no CI)", ylim = ylim_zoom, CI = FALSE, time.bins.line = TRUE)

plot.results.dispRity(out2$object, method = 2, colors = colors, colors_CI = colors_CI,
  main = "Duration", ylim = ylim, CI = TRUE)
plot.results.dispRity(out2$object, method = 2, colors = colors, colors_CI = colors_CI,
  main = "Duration (no CI)", ylim = ylim_zoom, CI = FALSE, time.bins.line = TRUE)

plot.results.dispRity(out2$object, method = 3, colors = colors, colors_CI = colors_CI,
  main = "Number", ylim = ylim, CI = TRUE)
plot.results.dispRity(out2$object, method = 3, colors = colors, colors_CI = colors_CI,
  main = "Number(no CI)", ylim = ylim_zoom, CI = FALSE, time.bins.line = TRUE)

Plotting the results histogram style

Another way for the plots would be kind of our histogram representation (you can certainly use ggplot here, below is just an example):

The results per method

## Set the list of available models
models <- list("equalbins", "acctran", "deltran", "punctuated", "gradual")

## Get the stratigraphic, duration and number data for the ages
strat_data_ages <- lapply(models, make.class.histogram, out1$results, method = 1)
durat_data_ages <- lapply(models, make.class.histogram, out1$results, method = 2)
numbe_data_ages <- lapply(models, make.class.histogram, out1$results, method = 3)

## Get the stratigraphic, duration and number data for the epochs
strat_data_epoc <- lapply(models, make.class.histogram, out2$results, method = 1)
durat_data_epoc <- lapply(models, make.class.histogram, out2$results, method = 2)
numbe_data_epoc <- lapply(models, make.class.histogram, out2$results, method = 3)

## The lower age and the lower disparity limit are arbitrarily fixed to 0 and 15 respectively.
lower_age <- 0; lower_lim <- 15
#TG: might want to change that for other datasets!

## limits functions
get.break <- function(x) return(x$breaks)
get.count <- function(x) return(x$counts)
## Ages limits
xlim_ages <- c(max(unlist(c(sapply(strat_data_ages, get.break), sapply(durat_data_ages, get.break),
  sapply(numbe_data_ages, get.break)))), lower_age)
ylim_ages <- c(lower_lim, max(unlist(c(sapply(strat_data_ages, get.count),
  sapply(durat_data_ages, get.count), sapply(numbe_data_ages, get.count)))))
## Epoch limits
xlim_epoc <- c(max(unlist(c(sapply(strat_data_epoc, get.break), sapply(durat_data_epoc, get.break),
  sapply(numbe_data_epoc, get.break)))), lower_age)
ylim_epoc <- c(lower_lim, max(unlist(c(sapply(strat_data_epoc, get.count),
  sapply(durat_data_epoc, get.count), sapply(numbe_data_epoc, get.count)))))

## Colors
colors <- palette()[1:5]
colors[1] <- "grey"

## Graphics
op <- par(bty = "n", mfrow = c(3,2))

## Plot the stratigraphy - ages
plot(strat_data_ages[[1]], xlim = xlim_ages, ylim = ylim_ages , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Stratigraphy - ages")
## Plot the slicing
for(slice in 2:5) {
  lines(strat_data_ages[[slice]]$breaks, strat_data_ages[[slice]]$counts, col = colors[slice])
}

## Adding the legend
legend(x = "bottomleft", legend = c("time bins", "acctran", "deltran", "punctuated", "gradual"),
  col = colors[1:5], lty = 1, bg = "white")

## Plot the stratigraphy - epoch
plot(strat_data_epoc[[1]], xlim = xlim_epoc, ylim = ylim_epoc , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Stratigraphy - Epoch")
## Plot the slicing
for(slice in 2:5) {
  lines(strat_data_epoc[[slice]]$breaks, strat_data_epoc[[slice]]$counts, col = colors[slice])
}


## Plot the duration - ages
plot(durat_data_ages[[1]], xlim = xlim_ages, ylim = ylim_ages , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Duration - ages")
## Plot the slicing
for(slice in 2:5) {
  lines(durat_data_ages[[slice]]$breaks, durat_data_ages[[slice]]$counts, col = colors[slice])
}

## Plot the duration - epoch
plot(durat_data_epoc[[1]], xlim = xlim_epoc, ylim = ylim_epoc , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Duration - Epoch")
## Plot the slicing
for(slice in 2:5) {
  lines(durat_data_epoc[[slice]]$breaks, durat_data_epoc[[slice]]$counts, col = colors[slice])
}


## Plot the number - ages
plot(numbe_data_ages[[1]], xlim = xlim_ages, ylim = ylim_ages , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Number - ages")
## Plot the slicing
for(slice in 2:5) {
  lines(numbe_data_ages[[slice]]$breaks, numbe_data_ages[[slice]]$counts, col = colors[slice])
}

## Plot the number - epoch
plot(numbe_data_epoc[[1]], xlim = xlim_epoc, ylim = ylim_epoc , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Number - Epoch")
## Plot the slicing
for(slice in 2:5) {
  lines(numbe_data_epoc[[slice]]$breaks, numbe_data_epoc[[slice]]$counts, col = colors[slice])
}

Some of the plots are a bit off set and everything but it gives us an idea

The results per model

Something like that but with transparency:

plot(strat_data_epoc[[1]], xlim = xlim_epoc, ylim = ylim_epoc , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[1], main = "Epoch")
plot(durat_data_epoc[[1]], xlim = xlim_epoc, ylim = ylim_epoc , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[2], main = "", add = TRUE)
plot(numbe_data_epoc[[1]], xlim = xlim_epoc, ylim = ylim_epoc , xlab = "Time (Mya)",
  ylab = "Sum of variance", col = colors[3], main = "", add = TRUE)